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ABSTRACT 


A numerical method is outlined for solving the equation which de- 
scribes the solar modulation of cosmic rays in models where interplane- 
tary conditions can vary with heliocentric latitude. As an illustration 
of the use of this method, it is shown how variations in the modulation 
with latitude could produce the small radial gradients in the intensity 
that have been observed from the Pioneers 10 and 11 spacecraft. 
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INTRODUCTION 


In Fisk (1971), a numerical method is outlined for solving the 
equation which governs the solar modulation of galactic cosmic rays in a 
spherically-symmetric model for the interplanetary medium. In the last 
few years this method has been used in numerous studies of the mod- 
ulation problem. It is the purpose of the present paper to extend this 
method to solve the modulation equation in models where interplanetary 
conditions vary with heliocentric latitude. 

The cosmic-ray flux that is observed in the solar equatorial plane 
may well depend sensitively on interplanetary conditions at other helio- 
centric latitudes. For example, observed cosmic rays may have entered 
the inner solar system over the solar poles, and then diffused across 
the mean magnetic field into the equatorial plane. The distance that 
must be travelled by particles following field lines from the inter- 
stellar medium to the inner solar system is much shorter over the poles. 
The magnetic field in this region should lie nearly in the heliocentric 
radial direction, whereas n the equatorial plane it is wound in a tight 
spiral. it is also conceivable that cosmic rays may tend to diffuse out 
of the region near the equatorial plane into the mid-latitude regions, 
above ana below the plane. Solar activity appears to be enhanced in the 
mid-latitude regions, and thus the cosmic-ray intensity there may be 
depressed . 

As an illustration of the use of the numerical method outlined 
here, it is shown how latitude variations in the modulation can affect 
the radial gradient in the intensity seen in the equatorial plane. In 
particular, it is shown how these variations could produce the small gra- 
dients which are measured from Pioneers 10 and 11 (Teegarden et a 1 . , 

1971; McKibben et al . , 1971; Van Allen, 1972 a, b). 
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THE METHOD 


Consider a model in which interplanetary conditions vary with 
heliocentric radial distance r and polar angle 0; the latter is measured 
relative to the axis of rotation of the sun. In steady-state conditions, 
the cosmic-ray omni-directional distribution function f (number of par- 
ticles per unit volume of phase space, averaged over particle direction) 
behaves in such a model according to the equation (Parker, 1965; Cleeson 
and Axford, 1967; Jokipii and Parker, 1970; Fisk a_l. , 1973): 

-£_ 2 3_(r 2 V)3f = l 2 i_(r 2 < 3f) + 1 3 (sin 9 k £) _ V 3f (1) 

3r 3r 3p r 3r r 3r 2 . ' 36 3r 

r sin9 39 

The distribution function f is related to the differential intensity j, 

2 

per unit interval of kinetic energy T, by j = f/p (Forman, 197C). 

The term on the left side of (1) and the third term on the right 
si<ie describe the effects of convection and adiabatic deceleration in 
the expanding solar wind. Here p is the magnitude of particle momentum, 
and V(r,9) is the solar wind speed. It is assumed that the solar wind 

flows only in the radial direction. The first and second terms on the 

right side of (1) describe the effects of diffusion in the radial and 

polar directions, respectively. The diffusion coefficient can also 

be expressed as 
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Here, tc„ is the diffusion coefficient for propagation parallel to the 
mean magnetic field, which is assumed to lie on cones of constant 6 and 
to make an angle ij> with the radial direction. The diffusion coefficient 
k ( describes the propagation normal to the mean field, along cones of 


2 


constant 0. In most cases, r equals the cross-field diffusion 

1 u 

coefficient for propagation in the polar direction. Gradient and cur- 
vature drifts of the particles in the large-scale interplanetary field 
are not considered here. However, the effects of such motion could be 
included in (1) simply by adding a term v^ * Vf, where v^ is the drift 
velocity. 

Equation (1) is a parabolic partial differential equation. As 
such, it can be readily solved by using the Crar.k-Nicholson implicit 
finite difference technique, which was developed for solving space-time 
diffusion equations. Momentum in (1) is the analogue of time in the 
simple diffusion equation. In spherically-symmetric models, where the 
second term on the right side of (1) is ignored, (1) can be solved by a 
straight-forward application of the Crank-Nicholson technique, as is 
done in Fisk (1971). In models where latitude variations are included, 
the appropriate technique is the alternating - direction modification of 
the Crank-Nicholson method, which was suggested by Douglas (1962) and 
elaborated on by Douglas and Gunn (1964) . A useful discussion of this 
method can be found in Carnahan et al . (1969) 

In finite-difference schemes a value for f is obtained at a series 
of grid points f (I, J, K) . Here 9 = I‘A0, r = J‘Ar + r^, and p = K‘Ap 
+ p^. The spacings between grid points are A9, Ar, and Ap; A0 = r/M, 
where M is the largest value of I. Boundary conditions must be speci- 
fied at I = 0 and I = M, corresponding to 0 = 0 and tt , t espectively . 
Similarly, bourdary conditions must be specified at r = r , the smal- 
lest radial distance considered, and at R =N‘Ar + r Q , the largest value 
of r. An "initial." condition must be given at p = p^, i.e. f (I, J, 0) 
must be given at all I and J. 
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In the simple Crank-Nicholson technique, derivatives with respect 
to 0 are replaced with finite difference equations which relate f (1-1, 

J, K), f (I, J, K) and f (I + 1, J, K) to f (1-1, J, K + 1) , f (I, J, 

K + 1) and f (I + 1, J, K + 1). Similarly, derivatives with respect to 
r relate f (I, J-l, K) , f (I, J, K) and f (I, J + 1, K) to their counter- 
parts at K + 1. Derivatives with respect to p involve f (I, J, K) and 

f (I, J, K + 1). By starting with the known values of f at h = 0, the 

value of f at larger K would then have to be determined in this scheme 

by solving at each step M*N simultaneous linear equations. For most 

applications this procedure is impractical. 

In the alternating-direction modification to the Crank-Nicholson 
techni 'ue, a two-step procedure is followed. In the first step the 
derivatives with respect to 0 are defined in terms of f (1-1, J, K) , 
f (I, J, K) and f (I + 1, J, K) , and intermediate values of f* (1-1, J, K 
+ 1) , f* (I, J, K + 1), and f* (I + 1, J, K + 1). The derivatives with 
respect to p also relate f (I, J, K) to f* (I, J, K + 1). The derivatives 
with respect to r, however, involve only f (I, J-l, K) , f (I, J, K) and 
f (I, J + 1, K). To determine f*, then, requires the solution of M 
simultaneous, linear equations, N times, which is a more tractable procedure. 
In the second step, the derivatives with respect to 3 remain defined in 

terms of f at K., and f*, but now the derivatives with respect to r and p 
are defined in terms of f at K and the required values of f at K + 1. 

"•o determine f at K + 1, N simultaneous, linear equations must be solved 
M times. By starting with f at K = 0, this two-step procedure is repeated 
at each value of K. A fast algorithm for solving the simultaneous 
equations can be found in Diaz (1958). 

For the boundary conditions at 0 = 0 and tt it is appropriate to 
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require simply that 3 f / 30 3 0. Some caution must be exercised in im- 
posing this condition, however, since (1) contains a term which is 
proportional to (l/tan0)3f/39. With an application of L'Hospital's 
rule, (l/tan0) 3 f /36 = 3 2 f / 30 ? at 6 = 0 and it. The maximum radial dis- 
tance R is chosen so that at r > R the modulation is negligibly small. 
The boundary condition at r = R is then f(I, N, K) = f^(p), where f Q (p) 
is the unmodulated, interstellar distribution function. Cosmic rays 
are assumed here to impinge on the solar cavity isotropically. For the 
bo ndary condition at r = r , which is assumed to be a value small 
cor pared with 1 AU, it is appropriate to require that 3f/9r = 0. 

T t should be noted that in Fisk (1971) the boundary condition at r 

= r was treated with considerable care. The differential number den- 
o 

1/2 1/2 

sity U was scaled by r or u = r U, and then u was set equal to 
zero at r = 0. In retrospect, this precaution was unnecessary. The 
solution near eaich, for example, is relatively insensitive to the boun- 
dary condition at r = r^. It suffices then simply to choose a condition 
which is easy to program, and which is not physically unrealistic. 

The initial momentum p is chosen so that at p>p the modulation is 

o o 

small. The initial condition is then f(I, J, 0) = f (p ). The solu- 

o o 

tion at values of p < p^ is determined by stepping down in momentum, 
i.e. Ap is taken to be negative. In this regard it is convenient to 
change variables in (1) from p to In p. The spacing between grid points 
is now Aln p, which is taken to be a constant. At large values of p, 
then, where there is little modulation, the steps in momentum ('vpAfcnp) 
are large. At small p, where f varies rapidly, the step size in p is 


automatically reduced. 


AN ILLUSTRATION EXAMPLE 


Measurements from the Pioneer 10 and 11 spacecraft have shown that 
the radial gradient of the cosmic-ray intensity is small. Teegarden 
et al . (1973) and McKibben et^ a£. (1973) find that the integral proton 
gradient (the gradient of protons with energies above ^ 60 MeV) is only 
about 4%/AU between 1 and 3 AU. Van Allen (1972 a, b) reports a gra- 
dient which is consistent with zero over this distance. McKibben et^ al. 
(1975) find that the integral gradient remains at 'v 4%/ AU out to at 
least 5 AU. McDonald et^ al . (1975) have argued recently that the cor- 
rections for Jovian electrons in their previous analysis, as well as in 
the analyses of the other experiments, were inadequate. With these 
corrections taken into proper account, McDonald e_t al. (1975) suggest 
that all previous gradient measurements could be reduced. 

At the radia. distances sampled by Pioneer the principle contri- 
bution to the integral proton intensity (protons above 'v 60 MeV in 
energy) comes from particles with energies ^ 1 GeV (e.g., McKibben et^ 
al . 1973). Equivalently, the observed integral gradient is a measure of 
the gradient of these high energy particles. In spherically-syrametric 
models for the interplanetary medium the high energy gradient can, in 
turn, be related to the radial diffusion coefficient by the relationship 

13,1 „ CV (3) 

J3r " < r 

where C = - 31n f/91n p is the Compton-Getting coefficient (Gleescn and 
Axford, 1968 a; Fisk and Axford, 1969, 1970; Fisk et^ al., 1973). 

The radial diffusion coefficient can, in principle, also be de- 
termined from the observed properties of the interplanetary magnetic 
field. In the last few years there has been contraversy as to exactly 
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how this calculation should bo performed (cf. Fisk et ah , 1974) . Most of 
the disagreement, however, is concerned with the behavior of particles 
at relatively low energies. For high energy particles, where the particle 
gyro-radius r^ is greater than, or on the order of the correlation 
length X of the field fluctuations, most theories appear to be in reason- 
able agreement (e.g. Jokipii, 1966; Klimas and Sandri, 1971). 

In the review article by Jokipii (1971), for example, the parallel 
diffusion coefficient is given as 


v r 2 B 2 

S 9_ 

P (k = 0) 


, for r > X 


(4) 


Here, P (k = 0) is the power density in fluctuations normal to the 
mean field direction, evaluated at zero wavenumber; v is particle speed; 
and is the mean field strength. The perpendicular diffusion coef- 
ficient in this limit is 
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for r > X 
g ^ 


(5) 


The observed power density at small wavenumbers yields P^ y (k = 9) 

2 2 

=1-2-10 gauss cm (See, for example, the power spectra given in 
Jokipii and Coleman (1968), or in Fisk and Sari (1973)). The mean field 
strength near r = 1 AU is B^ = A.5'10 ^ gauss. Thus, (4) and (5) yield 
respectively 



<„ = 2-10 21 BR 2 

(6) 

and 

k , = 1.2* 10 2 1 6 

(7) 


where R is particle rigidity in units of GV, and £ = v/c, with c the 
speed of light. These values for « „ and k ( agree closely with the 
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results given In Jokipii (1971). With a typical correlation length of 
\ *v 1.5*10^071 (Jokipii and Coleman 1968; Fisk and Sari, 1973), (6) and 
(7) should be valid for particles with rigidities R ^ 2 GV, or equiva- 
lently, for protons with energies T > 1 GeV. The solar wind speed is V 

- 400 km/sec. Thus, by use of (3), (6) and (7) predict a gradient for 

1-CeV protons, or equivalently an integral gradient, of v, 15%/AU. 

Conceivably, (6) and (7) could be small by a factor of 4. These 

expressions are evaluated here at T = 1 GeV, which is at the lower limit 

of their range of validity. An observed integral gradient of 'v 4%/AU, 
is then not necessarily inconsistent with theory. However, if the in- 
tegral gradient is in fact small compared with 4%/AU, as is indicated in 
the observations of Van Allen (1972 a, b) or is implied in the suggestion 
of McDonald £t a_l. (1975) , then there appears to be a conflict between 
the theory which is applied in sphericallv-symmetric models, and ob- 
servation . 

One possible solution to this problem is to disgard the assumption 

of a sphericallv-symmetric interplanetary medium. If there are more 

particles above and below the solar equatorial plane, and these particles 

tend to diffuse into the region of the plane preferentially at small 

radial distances, such a process reduces the radial gradient. The 

spiral pattern of the interplanetary field is, in fact, less tightly 

wound above and below the equatorial plane. Among other results, this 

2 

effect increases the term cos ^ in (2) at higher latitudes, and thus it 
may reduce the modulation here, and increase the particle density. The 
particles above and below the plane will diffuse onto the plane pre- 
ferentially at small r, provided that k. increases less rapidly than 

u 

2 

proportional to r , as can be seen from (1). 
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There i.3 another criterium which must be satisfied for a latitude- 


dependent modulation model to be acceptable. The observed proton spec- 
trum has a slope near unity at energies < 80 MeV, i.e. 31nj/3lnT 5 1 or 
'lnf/31n p=0 (e.g. Rygg and Earl, )971). In this energy range, then, 
the terms on the right of (1) must all sum to nearly zero, or they must 
all be individually small. For the case considered here, a balance of 
the terms docs not seem to be possible. Particles are diffusing into 
the region near the equatorial plane, and thus the second term on the 
right of (1) is positive. If, after entering the equatorial plane, 
these particles tend to diffuse outward in radial distance, the first 
term in (1) can be negative. However, since particles also diffuse 
inward along the equatorial plane from the interstellar medium, this 
first term can never cancel the second. Circumstances in which the 
first and third terms can be made individuall) small are discussed in 
detail in Fisk et^ al . (1973). The second term can bt made adequately 
small by requiring that becomes small a w low energies. 

Consider then the following model. For particles with r^ > A 
<„ and tc | , as are given in (4) and (Z), are used. For particles with 
< \ , ic„ is taken to be the form given in Jokipii (1971), which is ex- 
pressed here as 



( 8 ) 


P (k - 01 r 
xx \ g, 


It is assumed in (8) that the power spectrum of field fluctuations falls 
-3/2 

off as k for large k. It is also recognized that (8) yields a value for 
k„ which is probably too small at low energies, and thus it will tend to 

overestimate the radial gradient. For example, with P^(k = 0) = 

22 -5 11 

1.5*10 Gauss cm, = 4.5*10 Gauss, and \ = 1.5*10 cm, (8) predicts 
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'v 0.015 AU. Solar 


that the mean free path for a 10 MeV proton is X 


mf p 


flare observations, however, appear to indicate that X ^ 0.1 All at 

mrp 

these energies (Ma Sung et al . , 1975). 

Particles are assumed to diffuse across the mian field direction by 
scattering off of fluctuations which have scale-si/es comparable to 


As can be seen from the formula given in Jokipii (1971), tf , is then 

v P (k - 0) ft \ 3/2 


= K 0 « 


2B * 


X 


, for r < X 

R - 


(^ 


o ' ' 

This form has the required behavior that becomes small at low energies. 

Cross-field diffusion in which particles follow random-walking 
f J -*ld lines across the mean field direction is thus ignored here. This 
• . Lar procest , which is discussed by Jokipii and Parker (1969), supposedly 
yields a diffusion coefficient that is roughtly equal to (7) at all 
energies. In practice, however, as is Indicated by observations of solar 
particles, k at low energies is nowhere near this large (Krimigis et^ al • , 


1971). 

2 

The quantity (k “ 0)/B is taken to be a constant equal to 

7.4'10^cm, independent of r and 6. (Note that near earth P (k=0) = 

2 2 5 

1.5*10 Gauss“cm and B = 4.5*10 Gauss.) The cross-field diffusion 

o 

2 

coefficient given in (5) is then constant, or equivalently r^/r varies 
2 

as 1/r . Higher anergy particles will thus enter the region of the 

equatorial plane preferentially at small r. This behavior for 
2 

P (k = 0)/B can result provided that the correlation length, and the 

XX o 

amplitude of the fluctuations relative to the mean field are both con- 
stants. Thus, X is taken to be 1.5*10 ^cm, independent of r and 0. The 
solar wind speed is taken to be similarly constant at V = 400km/sec. 

The mean magnetic field is assumed to execute the Archimedes spiral 
pattern appropriate for each latitude (Parker, 1963). The field is 
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normalized so that at r * 1 AU, B = 4.5*10 ^ Gauss, and it = 45°. In 

o 

the model considered here the variation in the field magnitude and 
direction with polar angle 9 is solely responsible for the latitude 
dependence of the modulation. 

The outer edge of the modulating region is placed at R = 25 AU. 

2 2 —1 3 

The unmodulated in“erstellar spectrum is taken to be f^ « (T_ + c p ) / 

p, where T^ is particle rest energy. With this form, the spectrum for 
the differential number density is given by a power law in total energy, 
with spectral index -2.6. The unmodulated spectrum is normalized to be 
equal to the unmodulated spectrum that was used by Urch an! Gleeson 
(1972) in their extensive study of the modulation. 

With the above parameters and upon assuming that the particles are 
protons, (1) has been solved numerically by using the technique outlined 
in the previous section. Shown in Figure 1 is the calculated differential 
proton intensity near earth. As can be seen here, this spectrum pro- 
vides a good fit to the proton spectrum that is observed in _972, at the 
beginning of the Pioneer 10 mission. The calculated spectrum also 
provides a reasonable fit to the spectra observed by McKibben et^ al . 

(1975) at later dates, from Pioneers 10 and 11. Also shown in Figure 1 
is the calculated differential radial gradient, determined here between 
the distances 1 and 5 AU, in the equatorial plane. Rather than yielding 
a gradient v 15%/AU, (6) and (7) in this latitude-dependent model yield 
a gradient for 1-GeV protons of < 1%/AU. The integral proton gradient 
(protons above 60 MeV in energy) is in .act 'w 0.6%/AU. Clearly, 
variations in the modulation with heliocentric latitude could be respon- 
sible for the small gradients seen from Pioneers 10 and 11. 

There is some additional evidence which supports the explanation 
for the small gradients given here. Axford et aJL. (1975) report that 
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the azimuthal anisotropy seen from Pioneers 10 and 11 for particles with 
energies > 480 MeV/nucleon is smaller than expected. The observed 
anisotropy is £ . = 0.46 + .11%, whereas the theorv, as it is applied in 

sphericallv-symmetric models, can vield C . > 1%. Axford et al. (1975) 

<P 

suggest that an unusually large ratio of tc /<„ for particles with energies 

> 1 GeV/nucleon could be responsible for the reduced anisotropy. An 

% 

alternative explanation is that particles diffusing onto the equatorial 
plane from higher latitudes, as in the model discussed here, cause an 

outward radial streaming which in turn diminishes £.. 

9 

As can be seen in the discussion in Fisk (1974), in steady-state 
conditions 



( 10 ) 


where is the radial anisotropy. It is assumed in (10) that ‘' ( /k„<<1. 

_2 

In the model given here, for example, < 10 for p? i> s rbove 

1 GeV in energy, at r > 3 AU. Also the azimuth? I and pol c. ^dients 
are assumed to be small; the latter assumption is appropriate since in 
the model given here the modulation is symmetric about the equatorial 
plane . 

In spherically-symmetric models, ^ = 0 for particles with energies 
greater than several hundred MeV/nucleon (Gleeson and Axford, 1968 b; 
Fisk and Axford, 1969). At r = 3 AU, for example, tan ij; = 3. With C = 

1 and V = 400 km/sec, (10) then yields in such models = 1.4% for 1 

<P 

GeV protons. In contrast, in the latitude dependent model considered 
here is positive, as can be seen in Figure 1. At r = 3 AU in the 
equatorial plane, £ r shown here is 'v 0.3% for 1 GeV protons. With the 
same values for tis C and V, (10) then yields in this case E - 0.48%. 

9 

in good agreement with the observations of Axford ejt a_l. (1975). 
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The radial anisotropy predicted here would be difficult to de- 
tect from the Pioneer spacecraft. The spin axes of these spacecraft 
point in the radial direction. Also, a detailed examination of the 
numerical solution obtained here reveals that the radial anisotropy is 
less important near 1 AU. Thus, this anisotropy is not expected to in- 
fluence the diurnal anisotropy that is seen by neutron monitors. 

CONCLUDING REMARKS 

In addition to providing an explanation for the small gradients 
seen from Pioneers 10 and 11, the example presented here also illustrates 
how sensitive the cosmic-ray intensity observed near earth may be to 
latitude-variations in the modulation. Recall that the variation in 
modulation with latitude in this model is caused solely by the expected 
latitude variations in the direction and magnitude of the interplanetary 
magnetic field. This effect alone, however, introduces considerable 
polar gradients in the cosmic-ray intensity, and as a result, signifi- 
cantly alters the flux seen near earth. Shown in Figure 2 is the cal- 
culated intensity at 1 AU, plotted as a function of polar angle 9, for 

two different energies. The curves marked k = 0 are solutions to (1) 

0 

with the parameters given in the previous section, but now with polar 
diffusion ignored. The variation in the modulation conditions with 
latitude in this case causes the low-energy intensity to vary by nearly 
three orders of magnitude between the polar end equational regions. The 

curves marked < 4 0 are plots of the solution obtained in the 

0 

previous section with polar diffusion included. Note in particular the 

25 MeV curve. The polar diffusion coefficient is extremely small in 

18 2 -1 

this case (tc. v 7’10 cm sec ), and yet polar diffusion can alter the 

0 

intensity near earth hy a factor ^ 5. 
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It is reasonable to expect that the cosmic-ray modulation observed 
near earth will be understood only when information is available on 
modulating conditions at other heliocentric latitudes. 
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FIGURE CAPTIONS 


A plot vs. kinetic energy of the differential intensity 
near earth, the differential radial gradient in the 
equatorial plane between 1 and 5 AU, and the differential 
radial anisotropy in the equatorial plane at 3 AU. 

These curves are obtained from the numerical solution 
to (1) which is discussed in the text. The data 
points are the observed proton intensity during 
quiet-times in mid - 1972. These points, which were 
obtained from the GSFC experiment on IMP-5, were 
kindly provided by M. A. Van Holiebeke. Iso shown 
is the unmodulated interstellar spectrum wnich was 
used in the numerical example. 

A plot vs. polar angle 6 of the differential intensity 
at 1 AU, for two energies: T = 1 GeV and 25 Me'/. 

The curves marked k = 0 are solutions to (1), 

U 

obtained by using the parameters described in the 
text, but with polar diffusion ignored. The curves 
marked i 0 are the corresponding solutions 
with polar diffusion included. 
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